%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 5 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Note: due to the use of a finite number of simulations 
%%% to evaluate (k1,k2,k3), the final results are random 
%%% and will not correspond exactly to those in the paper.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Define variables
nobs = 10^3;%We use nobs=10^6 to create Figure 5 in the paper
gam = 1;
T = 100;
M = eye(T)-(ones(T,1)*ones(1,T))./T;
nuvec = 3:0.1:10;

%%% (\psi,\theta) = (0.2,0.3)
psi2 = 0.2^2;
theta2 = 0.3^2;
theta2g = theta2-psi2;
%N=10
N = 10;
EUloss2fN10theta03 = zeros(length(nuvec),1);
EUloss3fN10psi02theta03 = zeros(length(nuvec),1);
j = 1;
for nu = nuvec
    k1vec = zeros(nobs,1);
    k2vec = zeros(nobs,1);
    k3vec = zeros(nobs,1);
    for i = 1:nobs
         lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
         B = lambda*M*lambda;
         oneslambda = lambda*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         k1vec(i,1) = trace(Mat1);
         k2vec(i,1) = trace(Mat2);
         k3vec(i,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    k1 = a1*mean(k1vec);
    k2 = a2*mean(k2vec);
    k3 = a3*mean(k3vec); 
    a4 = (T-N-1)*(T-N-4)/(T*(T-2));
    %Two-fund rule
    cnorm = a4*theta2/(theta2+N/T);
    cellip = a4*k1*theta2/(k2*theta2+k3*N/T);
    EUnormal = (1/(2*gam))*(T/(T-N-2))*(2*cnorm*k1*theta2-...
               cnorm^2*(k2*theta2+k3*N/T)/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*cellip*k1*theta2-...
               cellip^2*(k2*theta2+k3*N/T)/a4);
    EUloss2fN10theta03(j,1) = EUnormal - EUellip;
    %Three-fund rule
    c1norm = a4*psi2/(psi2+N/T);
    c2norm = a4*(N/T)/(psi2+N/T); 
    c1ellip = a4*k1*psi2/(k2*psi2+k3*N/T);
    c2ellip = a4*(k3/k2)*(k1*N/T)/(k2*psi2+k3*N/T);
    EUnormal =  (1/(2*gam))*(T/(T-N-2))*(2*c1norm*k1*theta2+2*c2norm*k1*theta2g-...
                c1norm^2*(k2*theta2+k3*N/T)/a4-c2norm^2*k2*theta2g/a4-...
                2*c1norm*c2norm*k2*theta2g/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*c1ellip*k1*theta2+2*c2ellip*k1*theta2g-...
               c1ellip^2*(k2*theta2+k3*N/T)/a4-c2ellip^2*k2*theta2g/a4-...
               2*c1ellip*c2ellip*k2*theta2g/a4);
    EUloss3fN10psi02theta03(j,1) = EUnormal - EUellip;
    j = j+1;
end
%N=30
N = 30;
EUloss2fN30theta03 = zeros(length(nuvec),1);
EUloss3fN30psi02theta03 = zeros(length(nuvec),1);
j = 1;
for nu = nuvec
    k1vec = zeros(nobs,1);
    k2vec = zeros(nobs,1);
    k3vec = zeros(nobs,1);
    for i = 1:nobs
         lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
         B = lambda*M*lambda;
         oneslambda = lambda*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         k1vec(i,1) = trace(Mat1);
         k2vec(i,1) = trace(Mat2);
         k3vec(i,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    k1 = a1*mean(k1vec);
    k2 = a2*mean(k2vec);
    k3 = a3*mean(k3vec); 
    a4 = (T-N-1)*(T-N-4)/(T*(T-2));
    %Two-fund rule
    cnorm = a4*theta2/(theta2+N/T);
    cellip = a4*k1*theta2/(k2*theta2+k3*N/T);
    EUnormal = (1/(2*gam))*(T/(T-N-2))*(2*cnorm*k1*theta2-...
               cnorm^2*(k2*theta2+k3*N/T)/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*cellip*k1*theta2-...
               cellip^2*(k2*theta2+k3*N/T)/a4);
    EUloss2fN30theta03(j,1) = EUnormal - EUellip;
    %Three-fund rule
    c1norm = a4*psi2/(psi2+N/T);
    c2norm = a4*(N/T)/(psi2+N/T); 
    c1ellip = a4*k1*psi2/(k2*psi2+k3*N/T);
    c2ellip = a4*(k3/k2)*(k1*N/T)/(k2*psi2+k3*N/T);
    EUnormal =  (1/(2*gam))*(T/(T-N-2))*(2*c1norm*k1*theta2+2*c2norm*k1*theta2g-...
                c1norm^2*(k2*theta2+k3*N/T)/a4-c2norm^2*k2*theta2g/a4-...
                2*c1norm*c2norm*k2*theta2g/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*c1ellip*k1*theta2+2*c2ellip*k1*theta2g-...
               c1ellip^2*(k2*theta2+k3*N/T)/a4-c2ellip^2*k2*theta2g/a4-...
               2*c1ellip*c2ellip*k2*theta2g/a4);
    EUloss3fN30psi02theta03(j,1) = EUnormal - EUellip;
    j = j+1;
end
%N=50
N = 50;
EUloss2fN50theta03 = zeros(length(nuvec),1);
EUloss3fN50psi02theta03 = zeros(length(nuvec),1);
j = 1;
for nu = nuvec
    k1vec = zeros(nobs,1);
    k2vec = zeros(nobs,1);
    k3vec = zeros(nobs,1);
    for i = 1:nobs
         lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
         B = lambda*M*lambda;
         oneslambda = lambda*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         k1vec(i,1) = trace(Mat1);
         k2vec(i,1) = trace(Mat2);
         k3vec(i,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    k1 = a1*mean(k1vec);
    k2 = a2*mean(k2vec);
    k3 = a3*mean(k3vec); 
    a4 = (T-N-1)*(T-N-4)/(T*(T-2));
    %Two-fund rule
    cnorm = a4*theta2/(theta2+N/T);
    cellip = a4*k1*theta2/(k2*theta2+k3*N/T);
    EUnormal = (1/(2*gam))*(T/(T-N-2))*(2*cnorm*k1*theta2-...
               cnorm^2*(k2*theta2+k3*N/T)/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*cellip*k1*theta2-...
               cellip^2*(k2*theta2+k3*N/T)/a4);
    EUloss2fN50theta03(j,1) = EUnormal - EUellip;
    %Three-fund rule
    c1norm = a4*psi2/(psi2+N/T);
    c2norm = a4*(N/T)/(psi2+N/T); 
    c1ellip = a4*k1*psi2/(k2*psi2+k3*N/T);
    c2ellip = a4*(k3/k2)*(k1*N/T)/(k2*psi2+k3*N/T);
    EUnormal =  (1/(2*gam))*(T/(T-N-2))*(2*c1norm*k1*theta2+2*c2norm*k1*theta2g-...
                c1norm^2*(k2*theta2+k3*N/T)/a4-c2norm^2*k2*theta2g/a4-...
                2*c1norm*c2norm*k2*theta2g/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*c1ellip*k1*theta2+2*c2ellip*k1*theta2g-...
               c1ellip^2*(k2*theta2+k3*N/T)/a4-c2ellip^2*k2*theta2g/a4-...
               2*c1ellip*c2ellip*k2*theta2g/a4);
    EUloss3fN50psi02theta03(j,1) = EUnormal - EUellip;
    j = j+1;
end

%%% (\psi,\theta) = (0.4,0.5)
psi2 = 0.4^2;
theta2 = 0.5^2;
theta2g = theta2-psi2;
%N=10
N = 10;
EUloss2fN10theta05 = zeros(length(nuvec),1);
EUloss3fN10psi04theta05 = zeros(length(nuvec),1);
j = 1;
for nu = nuvec
    k1vec = zeros(nobs,1);
    k2vec = zeros(nobs,1);
    k3vec = zeros(nobs,1);
    for i = 1:nobs
         lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
         B = lambda*M*lambda;
         oneslambda = lambda*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         k1vec(i,1) = trace(Mat1);
         k2vec(i,1) = trace(Mat2);
         k3vec(i,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    k1 = a1*mean(k1vec);
    k2 = a2*mean(k2vec);
    k3 = a3*mean(k3vec); 
    a4 = (T-N-1)*(T-N-4)/(T*(T-2));
    %Two-fund rule
    cnorm = a4*theta2/(theta2+N/T);
    cellip = a4*k1*theta2/(k2*theta2+k3*N/T);
    EUnormal = (1/(2*gam))*(T/(T-N-2))*(2*cnorm*k1*theta2-...
               cnorm^2*(k2*theta2+k3*N/T)/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*cellip*k1*theta2-...
               cellip^2*(k2*theta2+k3*N/T)/a4);
    EUloss2fN10theta05(j,1) = EUnormal - EUellip;
    %Three-fund rule
    c1norm = a4*psi2/(psi2+N/T);
    c2norm = a4*(N/T)/(psi2+N/T); 
    c1ellip = a4*k1*psi2/(k2*psi2+k3*N/T);
    c2ellip = a4*(k3/k2)*(k1*N/T)/(k2*psi2+k3*N/T);
    EUnormal =  (1/(2*gam))*(T/(T-N-2))*(2*c1norm*k1*theta2+2*c2norm*k1*theta2g-...
                c1norm^2*(k2*theta2+k3*N/T)/a4-c2norm^2*k2*theta2g/a4-...
                2*c1norm*c2norm*k2*theta2g/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*c1ellip*k1*theta2+2*c2ellip*k1*theta2g-...
               c1ellip^2*(k2*theta2+k3*N/T)/a4-c2ellip^2*k2*theta2g/a4-...
               2*c1ellip*c2ellip*k2*theta2g/a4);
    EUloss3fN10psi04theta05(j,1) = EUnormal - EUellip;
    j = j+1;
end
%N=30
N = 30;
EUloss2fN30theta05 = zeros(length(nuvec),1);
EUloss3fN30psi04theta05 = zeros(length(nuvec),1);
j = 1;
for nu = nuvec
    k1vec = zeros(nobs,1);
    k2vec = zeros(nobs,1);
    k3vec = zeros(nobs,1);
    for i = 1:nobs
         lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
         B = lambda*M*lambda;
         oneslambda = lambda*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         k1vec(i,1) = trace(Mat1);
         k2vec(i,1) = trace(Mat2);
         k3vec(i,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    k1 = a1*mean(k1vec);
    k2 = a2*mean(k2vec);
    k3 = a3*mean(k3vec); 
    a4 = (T-N-1)*(T-N-4)/(T*(T-2));
    %Two-fund rule
    cnorm = a4*theta2/(theta2+N/T);
    cellip = a4*k1*theta2/(k2*theta2+k3*N/T);
    EUnormal = (1/(2*gam))*(T/(T-N-2))*(2*cnorm*k1*theta2-...
               cnorm^2*(k2*theta2+k3*N/T)/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*cellip*k1*theta2-...
               cellip^2*(k2*theta2+k3*N/T)/a4);
    EUloss2fN30theta05(j,1) = EUnormal - EUellip;
    %Three-fund rule
    c1norm = a4*psi2/(psi2+N/T);
    c2norm = a4*(N/T)/(psi2+N/T); 
    c1ellip = a4*k1*psi2/(k2*psi2+k3*N/T);
    c2ellip = a4*(k3/k2)*(k1*N/T)/(k2*psi2+k3*N/T);
    EUnormal =  (1/(2*gam))*(T/(T-N-2))*(2*c1norm*k1*theta2+2*c2norm*k1*theta2g-...
                c1norm^2*(k2*theta2+k3*N/T)/a4-c2norm^2*k2*theta2g/a4-...
                2*c1norm*c2norm*k2*theta2g/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*c1ellip*k1*theta2+2*c2ellip*k1*theta2g-...
               c1ellip^2*(k2*theta2+k3*N/T)/a4-c2ellip^2*k2*theta2g/a4-...
               2*c1ellip*c2ellip*k2*theta2g/a4);
    EUloss3fN30psi04theta05(j,1) = EUnormal - EUellip;
    j = j+1;
end
%N=50
N = 50;
EUloss2fN50theta05 = zeros(length(nuvec),1);
EUloss3fN50psi04theta05 = zeros(length(nuvec),1);
j = 1;
for nu = nuvec
    k1vec = zeros(nobs,1);
    k2vec = zeros(nobs,1);
    k3vec = zeros(nobs,1);
    for i = 1:nobs
         lambda = diag(sqrt((nu-2)./chi2rnd(nu,[T,1])));
         B = lambda*M*lambda;
         oneslambda = lambda*ones(T,1);
         Y = randn(T,N);
         Mat1 = inv(Y'*B*Y);
         Mat2 = Mat1^2;
         k1vec(i,1) = trace(Mat1);
         k2vec(i,1) = trace(Mat2);
         k3vec(i,1) = oneslambda'*Y*Mat2*Y'*oneslambda;
    end
    a1 = (T-N-2)/N;
    a2 = a1*(T-N-1)*(T-N-4)/(T-2);
    a3 = a2/T;
    k1 = a1*mean(k1vec);
    k2 = a2*mean(k2vec);
    k3 = a3*mean(k3vec); 
    a4 = (T-N-1)*(T-N-4)/(T*(T-2));
    %Two-fund rule
    cnorm = a4*theta2/(theta2+N/T);
    cellip = a4*k1*theta2/(k2*theta2+k3*N/T);
    EUnormal = (1/(2*gam))*(T/(T-N-2))*(2*cnorm*k1*theta2-...
               cnorm^2*(k2*theta2+k3*N/T)/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*cellip*k1*theta2-...
               cellip^2*(k2*theta2+k3*N/T)/a4);
    EUloss2fN50theta05(j,1) = EUnormal - EUellip;
    %Three-fund rule
    c1norm = a4*psi2/(psi2+N/T);
    c2norm = a4*(N/T)/(psi2+N/T); 
    c1ellip = a4*k1*psi2/(k2*psi2+k3*N/T);
    c2ellip = a4*(k3/k2)*(k1*N/T)/(k2*psi2+k3*N/T);
    EUnormal =  (1/(2*gam))*(T/(T-N-2))*(2*c1norm*k1*theta2+2*c2norm*k1*theta2g-...
                c1norm^2*(k2*theta2+k3*N/T)/a4-c2norm^2*k2*theta2g/a4-...
                2*c1norm*c2norm*k2*theta2g/a4);
    EUellip =  (1/(2*gam))*(T/(T-N-2))*(2*c1ellip*k1*theta2+2*c2ellip*k1*theta2g-...
               c1ellip^2*(k2*theta2+k3*N/T)/a4-c2ellip^2*k2*theta2g/a4-...
               2*c1ellip*c2ellip*k2*theta2g/a4);
    EUloss3fN50psi04theta05(j,1) = EUnormal - EUellip;
    j = j+1;
end

%%% Create Figure 5
fig = figure;
%Subplot 1
subplot(3,2,1)
plot(nuvec,12*EUloss2fN10theta03,'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(nuvec,12*EUloss2fN30theta03,'--b','LineWidth',2);
plot(nuvec,12*EUloss2fN50theta03,':r','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
legend('$N=10$','$N=30$','$N=50$','Location','southeast','interpreter','latex');
grid on
title('Two-fund rule, $\theta=0.3$','interpreter','latex');
xlabel('$\nu$','interpreter','latex');
xlim([3 10]);
ylim([-0.03 0]);
xticks([3 4 5 6 7 8 9 10]);
yticks([-0.03 -0.02 -0.01 0]);
%Subplot 2
subplot(3,2,2)
plot(nuvec,12*EUloss3fN10psi02theta03,'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(nuvec,12*EUloss3fN30psi02theta03,'--b','LineWidth',2);
plot(nuvec,12*EUloss3fN50psi02theta03,':r','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('Three-fund rule, $(\psi,\theta)=(0.2,0.3)$','interpreter','latex');
xlabel('$\nu$','interpreter','latex');
xlim([3 10]);
ylim([-0.3 0]);
xticks([3 4 5 6 7 8 9 10]);
yticks([-0.3 -0.2 -0.1 0]);
%Subplot 3
subplot(3,2,3)
plot(nuvec,12*EUloss2fN10theta05,'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(nuvec,12*EUloss2fN30theta05,'--b','LineWidth',2);
plot(nuvec,12*EUloss2fN50theta05,':r','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('Two-fund rule, $\theta=0.5$','interpreter','latex');
xlabel('$\nu$','interpreter','latex');
xlim([3 10]);
ylim([-0.3 0]);
xticks([3 4 5 6 7 8 9 10]);
yticks([-0.3 -0.2 -0.1 0]);
%Subplot 4
subplot(3,2,4)
plot(nuvec,12*EUloss3fN10psi04theta05,'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(nuvec,12*EUloss3fN30psi04theta05,'--b','LineWidth',2);
plot(nuvec,12*EUloss3fN50psi04theta05,':r','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('Three-fund rule, $(\psi,\theta)=(0.4,0.5)$','interpreter','latex');
xlabel('$\nu$','interpreter','latex');
xlim([3 10]);
ylim([-0.6 0]);
xticks([3 4 5 6 7 8 9 10]);
yticks([-0.6 -0.4 -0.2 0]);
fontsize(fig,16,"points");